# Load libraries, install if needed
pkgs <- c("tidyverse", "tidylog", "RCurl", "devtools", "patchwork",
"viridis", "RColorBrewer", "here", "raster", "concaveman")
if(length(setdiff(pkgs,rownames(installed.packages()))) > 0){
install.packages(setdiff(pkgs, rownames(installed.packages())), dependencies = T)
}
invisible(lapply(pkgs, library, character.only = T))
# Set path
home <- here::here()
# Source code for map plots
# For some reason I can't source it directly from github...
# You need: # devtools::install_github("seananderson/ggsidekick") # not on CRAN; library(ggsidekick)
# devtools::source_url("https://raw.githubusercontent.com/maxlindmark/pred-prey-overlap/main/R/functions/map-plot.R")
source(paste0(home, "/R/functions/map-plot.R"))
options(ggplot2.continuous.colour = "viridis")Calculate climate trends
Intro
This scripts calculates trends in biomass and climate variables from climate-agnostic sdm-predictions (“04-fit-sdms-random.qmd”)
Load packages & source functions
Read data
# Read prediction grids
d <- read_csv(paste0(home, "/output/data_pred_grids_04_random_sdms.csv")) |>
separate(group, c("species", "life_stage"), sep = "_", remove = FALSE)Rows: 5660568 Columns: 42
── Column specification ────────────────────────────────────────────────────────
Delimiter: ","
chr (6): substrate, month_year, ices_rect, group, species, life_stage
dbl (36): X, Y, year, lon, lat, depth, quarter, month, oxy, temp, sal, sub_d...
ℹ Use `spec()` to retrieve the full column specification for this data.
ℹ Specify the column types or set `show_col_types = FALSE` to quiet this message.
Trim the prediction grid by species cumulative biomass
# Cumulative filter
# d_test <- data.frame(y = round(rnorm(10, 10, 2))) |>
# arrange(desc(y)) |>
# mutate(cumsum = cumsum(y),
# sum = sum(y),
# perc = cumsum / sum)
#
# d_test
#
# d_test |> filter(perc < 0.95)
d_grid_keep <- d |>
mutate(group = paste(group, quarter, sep = "_")) |>
group_by(X, Y, group) |>
# summarise across all years, but within group and quarter
summarise(grid_mean = mean(est)) |>
ungroup() |>
group_by(group) |>
arrange(desc(grid_mean)) |>
mutate(cumsum = cumsum(grid_mean),
sum = sum(grid_mean),
perc = cumsum / sum) |>
ungroup() |>
filter(perc < 0.95) |>
mutate(id = paste(X, Y, group, sep = "_"))mutate: changed 5,660,568 values (100%) of 'group' (0 new NA)
group_by: 3 grouping variables (X, Y, group)
summarise: now 195,192 rows and 4 columns, 2 group variables remaining (X, Y)
ungroup: no grouping variables
group_by: one grouping variable (group)
mutate (grouped): new variable 'cumsum' (double) with 195,192 unique values and 0% NA
new variable 'sum' (double) with 12 unique values and 0% NA
new variable 'perc' (double) with 195,181 unique values and 0% NA
ungroup: no grouping variables
filter: removed 21,973 rows (11%), 173,219 rows remaining
mutate: new variable 'id' (character) with 173,219 unique values and 0% NA
# Filter the main dataset to include these grid cells only
pred_grid_trim <- d |>
mutate(group = paste(group, quarter, sep = "_")) |>
mutate(id = paste(X, Y, group, sep = "_")) |>
filter(id %in% unique(d_grid_keep$id))mutate: changed 5,660,568 values (100%) of 'group' (0 new NA)
mutate: new variable 'id' (character) with 195,192 unique values and 0% NA
filter: removed 637,217 rows (11%), 5,023,351 rows remaining
# Filter id's that are in the trends data frame before summarizing... because that is also trimmed by trends-quantiles.
# So first trim by quantilesCalculate trends and trim quantiles.
# This is if I use a spatial trend model
biom_yr_slopes <- pred_grid_trim |>
mutate(id = paste(X, Y, group, sep = "_")) %>%
split(.$id) |>
purrr::map(~lm(exp(est) ~ year, data = .x)) |>
purrr::map_df(broom::tidy, .id = 'id') |>
filter(term == 'year') |>
dplyr::select(id, estimate) |>
rename(biom_slope = estimate) |>
separate(id, sep = "_", into = c("X", "Y", "species", "life_stage", "quarter"), remove = FALSE) |>
mutate(X = as.numeric(X),
Y = as.numeric(Y))mutate: no changes
filter: removed 173,219 rows (50%), 173,219 rows remaining
rename: renamed one variable (biom_slope)
mutate: converted 'X' from character to double (0 new NA)
converted 'Y' from character to double (0 new NA)
temp_yr_slopes <- pred_grid_trim |>
mutate(id = paste(X, Y, group, sep = "_")) %>%
split(.$id) |>
purrr::map(~lm(temp ~ year, data = .x)) |>
purrr::map_df(broom::tidy, .id = 'id') |>
filter(term == 'year') |>
dplyr::select(id, estimate) |>
rename(temp_slope = estimate) |>
separate(id, sep = "_", into = c("X", "Y", "species", "life_stage", "quarter"), remove = FALSE) |>
mutate(X = as.numeric(X),
Y = as.numeric(Y))mutate: no changes
filter: removed 173,219 rows (50%), 173,219 rows remaining
rename: renamed one variable (temp_slope)
mutate: converted 'X' from character to double (0 new NA)
converted 'Y' from character to double (0 new NA)
oxy_yr_slopes <- pred_grid_trim |>
mutate(id = paste(X, Y, group, sep = "_")) %>%
split(.$id) |>
purrr::map(~lm(oxy ~ year, data = .x)) |>
purrr::map_df(broom::tidy, .id = 'id') |>
filter(term == 'year') |>
dplyr::select(id, estimate) |>
rename(oxy_slope = estimate) |>
separate(id, sep = "_", into = c("X", "Y", "species", "life_stage", "quarter"), remove = FALSE) |>
mutate(X = as.numeric(X),
Y = as.numeric(Y))mutate: no changes
filter: removed 173,219 rows (50%), 173,219 rows remaining
rename: renamed one variable (oxy_slope)
mutate: converted 'X' from character to double (0 new NA)
converted 'Y' from character to double (0 new NA)
phi_yr_slopes <- pred_grid_trim |>
mutate(id = paste(X, Y, group, sep = "_")) %>%
split(.$id) |>
purrr::map(~lm(phi ~ year, data = .x)) |>
purrr::map_df(broom::tidy, .id = 'id') |>
filter(term == 'year') |>
dplyr::select(id, estimate) |>
rename(phi_slope = estimate) |>
separate(id, sep = "_", into = c("X", "Y", "species", "life_stage", "quarter"), remove = FALSE) |>
mutate(X = as.numeric(X),
Y = as.numeric(Y))mutate: no changes
filter: removed 173,219 rows (50%), 173,219 rows remaining
rename: renamed one variable (phi_slope)
mutate: converted 'X' from character to double (0 new NA)
converted 'Y' from character to double (0 new NA)
# Left join these data!
trends <- left_join(temp_yr_slopes, biom_yr_slopes, by = c("id", "X", "Y", "species", "life_stage", "quarter")) |>
left_join(oxy_yr_slopes, by = c("id", "X", "Y", "species", "life_stage", "quarter")) |>
left_join(phi_yr_slopes, by = c("id", "X", "Y", "species", "life_stage", "quarter"))left_join: added one column (biom_slope)
> rows only in x 0
> rows only in y ( 0)
> matched rows 173,219
> =========
> rows total 173,219
left_join: added one column (oxy_slope)
> rows only in x 0
> rows only in y ( 0)
> matched rows 173,219
> =========
> rows total 173,219
left_join: added one column (phi_slope)
> rows only in x 0
> rows only in y ( 0)
> matched rows 173,219
> =========
> rows total 173,219
# Trim slopes by quantiles
trends <- trends |>
group_by(species, life_stage, quarter) |>
mutate(upr_temp = quantile(temp_slope, probs = 0.995),
lwr_temp = quantile(temp_slope, probs = 0.005),
upr_oxy = quantile(oxy_slope, probs = 0.995),
lwr_oxy = quantile(oxy_slope, probs = 0.005),
upr_biom = quantile(biom_slope, probs = 0.995),
lwr_biom = quantile(biom_slope, probs = 0.005),
upr_phi = quantile(phi_slope, probs = 0.995),
lwr_phi = quantile(phi_slope, probs = 0.005)) |>
ungroup() |>
filter(temp_slope > lwr_temp & temp_slope < upr_temp) |>
filter(oxy_slope > lwr_oxy & oxy_slope < upr_oxy) |>
filter(biom_slope > lwr_biom & biom_slope < upr_biom) |>
filter(phi_slope > lwr_phi & phi_slope < upr_phi)group_by: 3 grouping variables (species, life_stage, quarter)
mutate (grouped): new variable 'upr_temp' (double) with 12 unique values and 0% NA
new variable 'lwr_temp' (double) with 12 unique values and 0% NA
new variable 'upr_oxy' (double) with 12 unique values and 0% NA
new variable 'lwr_oxy' (double) with 12 unique values and 0% NA
new variable 'upr_biom' (double) with 12 unique values and 0% NA
new variable 'lwr_biom' (double) with 12 unique values and 0% NA
new variable 'upr_phi' (double) with 12 unique values and 0% NA
new variable 'lwr_phi' (double) with 12 unique values and 0% NA
ungroup: no grouping variables
filter: removed 1,746 rows (1%), 171,473 rows remaining
filter: removed 1,746 rows (1%), 169,727 rows remaining
filter: removed 1,660 rows (1%), 168,067 rows remaining
filter: removed 245 rows (<1%), 167,822 rows remaining
Continue calculating average cell values
# filter the area-occupied-trimmed grid by the quantile-trimmed trends data frame
pred_grid_trim <- pred_grid_trim |>
filter(id %in% unique(trends$id))filter: removed 156,513 rows (3%), 4,866,838 rows remaining
# ... Now that we have the final set of grid cells by species (after trimming by core area occupied and slope-quantiles), summarize grid cell average oxygen, temperature, and biomass
pred_grid_sum <- pred_grid_trim |>
group_by(id) |>
summarise(mean_log_biomass = mean(est),
mean_temp = mean(temp),
mean_oxy = mean(oxy),
mean_phi = mean(phi)) |>
# Center the average across year variables
mutate(mean_log_biomass_ct = mean_log_biomass - mean(mean_log_biomass),
mean_temp_ct = mean_temp - mean(mean_temp),
mean_oxy_ct = mean_oxy - mean(mean_oxy),
mean_phi_ct = mean_phi - mean(mean_phi)) |>
ungroup() |>
dplyr::select(id, mean_log_biomass, mean_temp, mean_oxy, mean_phi,
mean_log_biomass_ct, mean_temp_ct, mean_oxy_ct, mean_phi_ct)group_by: one grouping variable (id)
summarise: now 167,822 rows and 5 columns, ungrouped
mutate: new variable 'mean_log_biomass_ct' (double) with 167,822 unique values and 0% NA
new variable 'mean_temp_ct' (double) with 31,997 unique values and 0% NA
new variable 'mean_oxy_ct' (double) with 31,996 unique values and 0% NA
new variable 'mean_phi_ct' (double) with 167,822 unique values and 0% NA
ungroup: no grouping variables
# Join in the average variables to the trends data
trends <- trends |>
mutate(id = paste(X, Y, species, life_stage, quarter, sep = "_")) |>
left_join(pred_grid_sum, by = "id")mutate: no changes
left_join: added 8 columns (mean_log_biomass, mean_temp, mean_oxy, mean_phi, mean_log_biomass_ct, …)
> rows only in x 0
> rows only in y ( 0)
> matched rows 167,822
> =========
> rows total 167,822
# Scale variables: which group?
trends <- trends |>
mutate(group = paste(species, life_stage, sep = "_")) |>
group_by(group) |>
mutate(temp_slope_sc = (temp_slope - mean(temp_slope)) / sd(temp_slope),
oxy_slope_sc = (oxy_slope - mean(oxy_slope)) / sd(oxy_slope),
phi_slope_sc = (phi_slope - mean(phi_slope)) / sd(phi_slope)) |>
ungroup() |>
mutate(quarter_f = as.factor(quarter))mutate: new variable 'group' (character) with 6 unique values and 0% NA
group_by: one grouping variable (group)
mutate (grouped): new variable 'temp_slope_sc' (double) with 167,822 unique values and 0% NA
new variable 'oxy_slope_sc' (double) with 167,822 unique values and 0% NA
new variable 'phi_slope_sc' (double) with 167,822 unique values and 0% NA
ungroup: no grouping variables
mutate: new variable 'quarter_f' (factor) with 2 unique values and 0% NA
# Fix names
trends <- trends |>
mutate(species = str_to_title(species),
life_stage = str_to_title(life_stage)) |>
mutate(group = paste(species, life_stage, sep = "_"))mutate: changed 167,822 values (100%) of 'species' (0 new NA)
changed 167,822 values (100%) of 'life_stage' (0 new NA)
mutate: changed 167,822 values (100%) of 'group' (0 new NA)
Plot trends!
# Because we trim the grids, let's add a polygon that is the overall outside border.
border_df <- trends
z <- concaveman(cbind(border_df$X, border_df$Y), concavity = 2)
plot(z)z_poly <- sp::SpatialPolygons(list(sp::Polygons(list(sp::Polygon(z)), ID = 1)))
domain <- as.data.frame(z_poly@polygons[[1]]@Polygons[[1]]@coords)
colnames(domain) <- c("X", "Y")
ggplot(domain, aes(x = X, y = Y)) +
geom_polygon(fill = NA, color = "tomato")# Plot velocities!
long_trends <- trends |>
# Not using phi here!
pivot_longer(c("biom_slope", "temp_slope", "oxy_slope")) |>
#pivot_longer(c("biom_slope", "temp_slope", "oxy_slope", "phi_slope")) |>
mutate(name2 = ifelse(name == "biom_slope", "Biotic", NA),
name2 = ifelse(name == "temp_slope", "Temperature", name2),
name2 = ifelse(name == "oxy_slope", "Oxygen", name2)#,
#name2 = ifelse(name == "phi_slope", "Metabolic index", name2)
) |>
mutate(loop_id = paste(paste0("Q", quarter), name2, sep = " ")) |>
mutate(unit = "(kg/decade)",
unit = ifelse(name2 == "Temperature", "(°C/decade)", unit),
unit = ifelse(name2 == "Oxygen", "(ml/L/decade)", unit)#,
#unit = ifelse(name2 == "Metabolic index", "(M.I./decade)", unit)
)
ggplot(long_trends, aes(value)) +
facet_wrap(~name, scales = "free") +
geom_histogram()# Loop through the variables... Note if oxygen, the tail is negative so ... trim by a small quantile
for(i in unique(long_trends$loop_id)) {
dd <- filter(long_trends, loop_id == i)
print(
plot_map_fc +
geom_raster(data = dd, aes(X*1000, Y*1000, fill = value)) +
facet_grid(life_stage ~ species) +
scale_fill_viridis(option = "rocket", name = paste(i, "trend", dd$unit[1])) +
geom_polygon(data = domain, aes(X*1000, Y*1000),
fill = NA, color = "red", size = 0.2)
)
# ggsave(paste0(home, "/figures/supp/trend_", str_replace(i, " ", "_"), ".pdf"), width = 17, height = 13, units = "cm")
}Save data
# Save data
write_csv(trends, paste0(home, "/data/clean/trends.csv"))